function [ VAR ] = irfvar( VAR ) 
%IRFVAR   IRFs for VAR under lower triangular Choleski ordering.
%
%  Inputs:
%    bet = reduced form VAR coefficients, dimension = # regressors x #variables
%  Sigma = VAR cov matrix of innovations
%    Perm= NyxNy permutation matrix to change causal ordering, i.e.
%          permutate rows of identity matrix.  Set to identity if causal
%          ordering is ordering of variables in y(t).
%    Ny  = # variables
%    p   = # lags
%   nirf = horizon for IRFs
% normalize = 1 if initial shock is normalized to unity, 0 if otherwise
%  Output:
%  virf  = nirf x Ny(#variables) x Ny (#shocks)       matrix

% Compute MA representation from companion form VAR:

Ny = VAR.Ny; p = VAR.p; P=VAR.Perm; Sigma=VAR.Sigma; bet = VAR.bet;
nirf = VAR.irhor; 

Sig_chol = P'*chol(P * Sigma * P')'*P ;                % transpose to get lower Cholesky factorization

B_comp = [ bet(1:p*Ny,:)' ; ...                       %companion form matrix of parameters 
          [eye(Ny*(p-1)),zeros(Ny*(p-1),Ny)]];
      
cc_dy   = eye(Ny*p);
msel    = [eye(Ny) ; zeros((p-1)*Ny,Ny)];

if VAR.normalize == 1;
    diagonal = diag(diag(Sig_chol));
    Sig_chol = Sig_chol*inv(diagonal);    % Hamilton p. 322: P = A D^(1/2) <=> A = P D^(-1/2), where Sigma = A D A', and A is a matrix with 1s in the principal diagonal
end

% note ordering: horizon, variable, shock
virf        = zeros(nirf,Ny,Ny); 
virf(1,:,:) = msel'*cc_dy*msel*Sig_chol;

for tt = 2:nirf;                 
    cc_dy = B_comp*cc_dy;
    virf(tt,:,:) = msel'*cc_dy*msel*Sig_chol;
end;
 
VAR.cholIRF = virf;
end

%B_comp looks like:
%   Phi_1   Phi_2 ... Phi_p
%    I       0          0
%    0       I     0    0
%    0       0     I    0


